#pragma once

#include "OrionFlames.h"
#include "GLVectorUtils.h"
#include "CoordinateMapping.h"

namespace OrionFlames
{
//A flame "Variation": a nonlinear function on 2D space. 
//Implemented both in C++ (mostly for testing/debugging)
//and in GLSL via the macros in MakeVariation.h.

typedef vec4 (*fn_taking_vec2_t)(const vec2 &t);

class OFAPI Variation
{
public:

	Variation() {}
	virtual ~Variation() { }
	
	//Map this location in 2D spac
	virtual vec2 apply(const vec2 &t)=0;
	
	//Return offical human-readable name of this Variation.
	virtual std::string getName() const =0;

	//Return GLSL-compatible code that reads a vec2 "t" 
	//and writes to the vec2 "p", applying this Variation.
	virtual std::string getCode()=0;
	
	//Call the function f on each of the inverse images 
	//of this point under this function.
	typedef vec4 (*fn_taking_vec2_t)(const vec2 &t);
	
	virtual vec4 sumInverse(const vec2 &p, fn_taking_vec2_t f)
		{return vec4(0.0);}
	
	//Return GLSL code that calls the function "f" with 
	//the vec2 "t" for each inverse of the vec2 "p", and
	//returns the vec4 sum of all the calls to "f".
	virtual std::string getInverseCode() {return "";}
	
	//Return the Jacobian determinant at this point.
	virtual float getArea(const vec2 &t)=0;
	
	//Return the GLSL code that returns the float jacobian
	//determinant, given the float source location (x,y).
	virtual std::string getAreaCode()=0;

protected:
	std::string m_VariationName;
#if 0 /* Do inverse and jacobian code verification at runtime */
* Silly inverse verification testing/printing code: */
Variation *curvar;
vec2 p;
vec4 f(const vec2 &t) {
	vec2 hp=curvar->apply(t);
	if (length(hp-p)>0.01) {
		printf("    bad %s inverse at point %f,%f  (hp=%f,%f, t=%f,%f)\n",
			curvar->getName().c_str(),
			p.x,p.y,
			hp.x,hp.y,
			t.x,t.y
		);
		//abort();
	}
	return vec4(1.001234);
}

* Verify that this Variation's inverse actually works: 
that I can start with random plane points, invert them,
and each of the inverted images maps back to the original point.
*/
void VerifyInverse(Variation *var,int seed) {
	curvar=var;
	Random48 r(seed);
	static int badcount=0;
	bool isbad=false;
	for (int rep=0;rep<10000;rep++) {
		p=vec2(20.0*r.nextFloat()-10.0,
			20.0*r.nextFloat()-10.0);
		if (var==&singleton_disc) /* *very* limited range for inverse */
			p=vec2(r.nextFloat()*0.7,r.nextFloat()*0.7);
		//Check inverse:
		vec4 v=var->sumInverse(p,f);
		
		//Check jacobian:
		float e=1.0e-4;
		vec2 t=p*0.1; /* limit range, to reduce roundoff */
		vec2 p0=var->apply(t);
		vec2 p1=var->apply(t+vec2(e,0));
		vec2 p2=var->apply(t+vec2(0,e));
		float area=fabs(p0.ccw(p1,p2)/(e*e));
		float jacobi=fabs(var->getArea(t));
		if (fabs((area-jacobi)/area)>0.1 &&
		    fabs(area-jacobi) > 1.0e-2) 
		{ /* Mismatch detected */
			printf("%d:    bad %s jacobian at point %f,%f  (area=%f, jacobi=%f)\n",
				badcount++,curvar->getName().c_str(),
				t.x,t.y,
				area,jacobi
			);
			isbad=true;
		}
		//else if (isbad) printf("."); /* indicating good jacobian */
	}
}

#else//no verification.
void VerifyInverse(int seed) { }//Ignore it.
#endif
};

//Linearly combine two nonlinear Variations.  Normally non-invertable. 
//This class would be useful for animations, but it's not used right now.
class Variationpair : public Variation
{
public:
	Variation *a,*b;
	float wa,wb;

	Variationpair(Variation *a_,float wa_,
		Variation *b_,float wb_)
		:a(a_),b(b_),wa(wa_),wb(wb_) {}
	virtual vec2 apply(const vec2 &src)
	{
		return wa*a->apply(src)+wb*b->apply(src);
	}

	virtual std::string getName() const
	{
		char nameStr[100];
		sprintf_s(nameStr,100,"(%f*%s+%f*%s)",
			wa,a->getName().c_str(),
			wb,b->getName().c_str());
		return nameStr;
	}

	virtual std::string getCode()
	{
		char weightStr[100];
		sprintf_s(weightStr,100,"%f*pa+%f*p\n;\n",wa,wb);
		return "vec2 pa;\n{" +
			a->getCode() +
			"} pa=p; {\n" +
			b->getCode() +
			"} p="+std::string(weightStr);
	}
};
}